Colloidal electrophoresis: Scaling analysis, 
Green— Kubo relation, and numerical results 

B Diinweg 1 , V Lobaskin 1,2 , K Seethalakshmy— Hariharan 1 and C 
Holm 1,3 

1 Max Planck Institute for Polymer Research, 
Ackermannweg 10, 55128 Mainz, Germany 

2 Physik-Department, Technische Universitat Miinchen, 
James-Franck-StraBe, 85747 Garching, Germany 

3 Frankfurt Institute for Advanced Studies (FIAS), J. W. Goethe-Universitat, 
Ruth-Moufang-Strafie 1, 60438 Frankfurt /Main, Germany 

E-mail: duenweg@mpip-mainz.mpg.de, lobaskin@tum.de, 
C . Holm@f ias . uni-f rankf urt . de 

Abstract. We consider electrophoresis of a single charged colloidal particle in a finite 
box with periodic boundary conditions, where added counterions and salt ions ensure 
charge neutrality. A systematic rescaling of the electrokinetic equations allows us 
to identify a minimum set of suitable dimcnsionless parameters, which, within this 
theoretical framework, determine the reduced electrophoretic mobility. It turns out 
that the salt-free case can, on the Mean Field level, be described in terms of just three 
parameters. A fourth parameter, which had previously been identified on the basis 
of straightforward dimensional analysis, can only be important beyond Mean Field. 
More complicated behavior is expected to arise when further ionic species are added. 
However, for a certain parameter regime, we can demonstrate that the salt-free case can 
be mapped onto a corresponding system containing additional salt. The Green-Kubo 
formula for the electrophoretic mobility is derived, and its usefulness demonstrated by 
simulation data. Finally, we report on finite-element solutions of the electrokinetic 
equations, using the commercial software package COMSOL. 
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1. Introduction 



The behavior of charge-stabilized colloidal dispersions in external electric fields is a 
classical topic of colloid physics [1]. Quantitative theoretical understanding is still 
incomplete today, although substantial progress has been achieved over the decades 
[2-9]. The difficulty lies in the complicated many-body nature of the problem, and 
hence only limiting cases are well understood. Beyond the physics of the "standard 
electrokinetic model" [6], which is essentially just a single-particle Mean Field theory 
(see below) , which nevertheless does describe a quite broad range of phenomena, current 
research focuses mainly on situations where this model is not applicable, or at least its 
applicability is not obvious. These include cases where non-Mean Field effects are 
important, i. e. higher valency or non-negligible size of the ions [7, 8, 10], or where 
the single-colloid picture is expected to break down, due to overlapping ionic clouds (or 
insufficient amount of screening by salt) [9]. This latter issue has also triggered detailed 
experiments [11-15] which measured the electrophoretic mobility in the low-salt regime. 
Furthermore, the problem has recently been studied by computer simulations [16-21]. 
The investigations of Ref. [21] were specifically targeted at the low-salt limit. The 
purpose of the present paper is to provide some more detailed theoretical and numerical 
background material which had to be omitted in Ref. [21]. We will start in Sec. 2 with 
a brief review of the simplest limiting cases of electrophoresis, followed by a summary 
of the observations made in Ref. [21]. The new material is then found in Sees. 3-5 
(briefly outlined at the end of Sec. 2), while we conclude in Sec. 6. 



2. Background: Review of limiting cases, and previous simulation results 

2.1. Single colloid without salt 

The simplest case of colloidal electrophoresis is obviously a single charged sphere of 
radius R and charge Ze (e denotes the positive elementary charge, and we assume 
Z > 0), immersed in a solvent of viscosity i] and dielectric constant e. Assuming zero 
salt concentration, and zero colloidal volume fraction $, the drift velocity v which results 
as a response to an external electric field E is just given by the Stokes formula, 

6nr]Rv = ZeE. (1) 

This is the so-called Hiickel limit [1, 3]. The electrophoretic mobility fi of the colloidal 
sphere, defined via 

v = fiE, (2) 

is hence given by 

We now introduce the zeta potential as the electrostatic potential at the colloid surface 
(with the understanding that it vanishes infinitely far away from the colloid): 
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This allows us to re-write (3) as 

^ = \k. (5) 

3 rj 

Based upon the thermal energy k B T as typical energy scale (k B denotes Boltzmann's 
constant and T the absolute temperature), we can introduce the dimensionless (reduced) 
zeta potential as 

(red = (6) 

On the other hand, the thermal energy, combined with electrostatics, provides a typical 
length scale, the Bjerrum length 

e 2 

which is nothing but the distance between two elementary charges such that their 
electrostatic energy is just k B T. This can be combined with the Stokes formula to 
define a useful mobility scale for electrokinetic phenomena: 

Defining the reduced mobility as 

l^red j (9) 

one sees that in the Hiickel limit one simply has 

Hred = (red- (10) 

2.2. Zeta potential vs. reduced charge 

In a more general context, the electrostatic potential at the surface of the colloid will of 
course no longer be given by (4). It will rather be diminished, as a result of the influence 
of the other charges. In order to clearly distinguish between the concepts of charge and 
potential, we will call 

~Z=^-^- = Z 1 ^ (11) 
4neRk B T R v ; 

the reduced (re-parametrized) charge (regardless of the physical situation), while the 
symbols ( and C, re d are reserved for the actual value of the surface electrostatic potential 
and its dimensionless counterpart. In the Hiickel limit (and only in this limit), ( re d and 
Z coincide. 



2.3. Screening 

An important aspect of electrophoresis is the screening of not only electrostatic, but also 
hydrodynamic interactions. As soon as one considers a system at a finite concentration, 
one has to take into account that it must be charge-neutral, at least on sufficiently 
large length scales: The charges (colloid charges plus ion charges) that are contained in 
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a sub-volume of linear dimension substantially larger than the colloid-colloid correlation 
length must add up to zero. The same is true (with arbitrary precision) in a computer 
simulation if one considers the simulation box as a whole (independently of the value of 
any correlation length). 

Now, the basic mechanism of hydrodynamic screening is the fact that the external 
electric field generates electric currents in both directions, which in turn generate 
hydrodynamic flows in both directions. In leading order, however, these flows cancel, 
since the total net force acting on the system (or sub-volume) is exactly zero. As a 
result, the flow around a moving charged colloid will not decay like 1/r (which would 
hold in the case of sedimentation, where the system responds to a gravitational field 
and the net force does not vanish), but much faster, ~ 1/r 3 [22]. The consequence is, on 
the one hand, that finite size effects in computer simulations are much less severe than 
in similar studies of sedimentation [10], and, on the other hand, that a single-particle 
picture will apply whenever the electrostatic interactions are sufficiently screened, as a 
result of high salt concentration. Indeed, it is well-known that in the high-salt limit 
the screening of electrostatics [1] is governed by the Debye length l D = k^ 1 , where the 
screening parameter k is proportional to the square root of the salt concentration c s : 



To be precise, (12) assumes monovalent salt ions, and c s denotes the total number of 
salt ions per unit volume (such that the number of ion pairs per unit volume is given by 
c s /2). Now, under conditions where l D is substantially smaller than the typical colloidal 
interparticle separation, it is clear that most of the space between the colloids is charge- 
neutral. Consequently, these regions are also force-free. In other words, in these regions 
there is no net flow, and all the hydrodynamic shear gradients and viscous dissipation 
processes are confined to the Debye layer as well. In this situation, one obviously can 
treat the problem in terms of a single-particle picture. However, even the problem of 
a single sphere surrounded by a charge cloud, with boundary condition of vanishing 
electrostatic potential, and finite salt concentration, for r — > oo, can in general be 
solved only numerically. This is the so-called "standard electrokinetic model" [6]. The 
reason for the mathematical difficulties is the non-linearity of the underlying Poisson- 
Boltzmann equation, which determines the ionic cloud structure. 

2.4- Smoluchowski limit 

A simple analytic solution is however possible in the limit of very high salt concentration 
such that l D <C R. Here the geometry is essentially planar, and one obtains the so-called 
Smoluchowski limit [1]: 

n = -C; (13) 

rj 

however, here the zeta potential is tiny, and in terms of the reduced charge one has 



k 2 = AttIbc, 



(12) 



f^red 




(14) 
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In the limit of infinitely strong screening (k — > oo), the salt completely shields the 
electric field from the particle, and correspondingly the mobility tends to zero. Of 
course, this is only a mathematical limit, which can never be reached in practice, since 
at a critical salt concentration the system of small ions will crystallize. Beyond this 
critical concentration the liquid-state Smoluchowski formula cannot work. 

2. 5. Simulations of the low-salt case 

While the case of high salt concentration can thus be considered as reasonably well 
understood (albeit in general only within the framework of numerics), a completely 
different situation arises when there is only little salt in the solution, or even none 
at all. In this case the ionic clouds are mainly formed by the counterions, and these 
will in general overlap. All the standard screening concepts, which are based upon 
assuming a decay of the electrostatic potential and of the charge density, on a length scale 
smaller than the colloid-colloid separation, are no longer expected to work. Nevertheless, 
namely due to the weak screening, some simplifying assumptions can still be made for 
suspensions of strongly charged colloids. As the colloids in this regime strongly repel 
each other, they are usually well ordered so that their minimal separation amounts to 
the mean interparticle distance d ~ i?$ -1 / 3 . Thus, the screening at r < d will be 
exclusively due to counterions and the phenomena that happen on this length scale will 
be governed by the mean counterion concentration. These ideas proved to be useful for 
describing static structure and colloidal interactions at low salt [23]. 

We have studied this case by computer simulations. In essence, our method is 
Molecular Dynamics (MD) for the charged colloid, the explicit (counter or salt) ions, 
and the solvent. However, for computational efficiency the latter is replaced by a lattice 
Boltzmann (LB) system which is coupled to the particles by a Stokes friction coefficient. 
This method, which has been designed as an efficient and easy way of simulating systems 
with hydrodynamic interactions, has been described in Refs. [24, 25], and is discussed in 
detail in a forthcoming review article [26]. Langevin noise is added to both the particles 
and the LB system to keep the temperature constant. The colloid is modeled as a 
"raspberry" [19, 20], i. e. a large central particle with a wrapping consisting of a tether 
of small particles. The most important results of this study have been communicated 
in Ref. [21], and can be summarized as follows: 

(i) fi re d is a dimensionless quantity, and hence can only depend on dimensionless 
parameters of the system. As a starting point, we have made no further theoretical 
assumptions. In the salt-free case, we can then identify four such parameters pi, ■ ■ ■ ,Pa, 
which we choose in such a way that two of them resemble most closely those quantities 
which have proven useful in the "salty" case: these are p\ = kR and p 2 = Z (cf. (10) 
and (14)). In the present case, however, k is not calculated from the salt concentration, 
but rather from the counterion concentration: 

k 2 = 4irl B c, (15) 
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with 



NZ 



Z 



3 



(16) 



c = 



4vri? 3 



where V is the system volume, and N the number of simulated colloids. Obviously, (15) 
and (16) imply the relation 



in other words, kR is nothing but a re-parametrized volume fraction. It should be 
emphasized that due to assumed strong charge asymmetry between the colloids and 
the counterions, which both constitute the screening medium, the resulting charge 
distribution is strongly inhomogeneous and the standard Debye screening concept cannot 
be implied. The remaining two scaling variables are p 3 = Ib/cl and p± = Ib/R, where a 
is the counterion radius. 

(ii) For a reasonable choice of parameters (Ib/R not too large, and l B /a of order 
unity, as is typical for ions in water) the dependence on p 3 and p<± can be ignored. 

(iii) In terms of Z and k,R, quite good agreement can be achieved with experiments, 
provided Z is replaced by an effective charge, calculated from charge renormalization 



(iv) Finite-size effects are weak, and hence one can obtain the data for a certain 
finite volume fraction by just simulating a single sphere in a suitably chosen finite box. 

(v) The effect of added salt is similar to that of increased volume fraction. It turns 
out that it is possible, within good approximation, to just combine these two effects into 
one single parameter 



which has a certain justification within a simplified linearized Poisson-Boltzmann 
theory [29]. 

The purpose of the present paper is to provide a theoretical background for the 
observations reported in Ref. [21] and derive some essential relations needed for the 
further analysis of the electrophoresis at finite colloidal concentrations. In particular, 
we describe our rescaling procedure in more detail. We feel that this will become 
particularly transparent when done in terms of the electrokinetic equations [1], which 
can be viewed as the Mean Field description of the system we have simulated — in 
contrast to the simulation, the counterions are not considered as discrete particles, but 
rather as concentration fields. Section 3 thus presents these equations, and outlines the 
rescaling procedure. An important result of this analysis is that the dependence on 
p 4 can indeed be ignored on the level of the electrokinetic equations — this parameter 
therefore describes deviations from Mean Field behavior, if there are any. Section 4 
discusses the problem of linear response, i. e. how to check that the non-equilibrium 
simulations employ a sufficiently weak external field. We have solved this by comparing 
the results with control calculations in strict thermal equilibrium, where the mobility 
was calculated by Green-Kubo integration. As far as we know, this formula has not yet 
been presented in the literature, and we will derive it here. Finally, in Sec. 5 we present 




2 



(17) 



[23, 27, 28]. 



k 2 = A-kIb (c + c s ) , 



(18) 
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some data which we have obtained by direct numerical solution of the elect rokinetic 
equations, using the commercial finite element package COMSOL 3.3. 

3. Rescaling of the electrokinetic equations 

In the stationary state, the electrokinetic equations are given by 

V-u = 0, (19) 
-Vp + r]V 2 v-e(y^)^2ziCi = 0, (20) 

i 

V • (-AVq - ^ e ^(W)q + vc)j = 0, (21) 
V 2 ^ + -eJ2 z i°i = °- ( 22 ) 

i 

Equation (19) is the incompressibility condition for the velocity field v, while (20) is 
the Stokes equation for zero Reynolds number flow, where the forces resulting from 
the hydrostatic pressure p and the viscous dissipation are balanced against the electric 
force. Here, \1/ denotes the electrostatic potential, while q is the concentration (number 
of particles per unit volume) of the ith ionic species. We will adopt the convention 
that i = corresponds to the counterions, while i > 1 denotes various types of salt 
ions. Each ion of species i carries a charge Zit. Hence, the total charge density is 
given by e^^Q; this term appears also in the Poisson equation for the electrostatic 
potential, (22), where the boundary conditions for ^ implicitly define the external 
driving field. Finally, (21) is the so-called Nernst-Planck equation (convection-diffusion 
equation) which describes the mass conservation of ionic species i. Here Di denotes the 
collective diffusion coefficient of species i, while Di/(kBT) is (via the Einstein relation) 
the corresponding ionic mobility. The ionic current consists of three contributions: the 
diffusion current, the drift relative to the surrounding solvent, induced by the electric 
force density —eziCi'V^ , and finally the convective current induced by the motion of the 
fluid. 

We now introduce 

Mi = J d 3 ra, (23) 

the number of ions of species % in the solution, where the integration extends over the 
finite volume V of the system. Obviously, the counterions just compensate the colloid 
charge, and hence we have 

-z M = Z; (24) 

note that Zq < 0, and we consider only a single colloid in the volume. Similarly, the 
charges of the salt ions compensate each other, and hence 

J>Mi = 0. (25) 

i>i 
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In the case without external driving, we have {7=0, and the Stokes equation reduces 
to an equation which determines the pressure. The Nernst-Planck equation, together 
with the Poisson equation, then just becomes the Poisson-Boltzmann equation: 

VhiQ + = 0, (26) 

V 2 ^ + 47^5^ = 0, (27) 

i 

where we have introduced the abbreviation 

* = i^r < 28 > 

In accordance with Ref. [29] and (18), we introduce the parameter 

V z 2 M- 

= 4^^A (29) 

where however no direct connection to a linearized Poisson-Boltzmann equation is 
implied. We now use R' 1 as our elementary unit of length and write 

V = reV, (30) 

which transforms the Poisson equation into a fully non-dimensional form: 

V 2 ^ + ^^q = 0, (31) 



where non-dimensional concentrations q are introduced via 

h. (32) 



Anl B V 

In these scaled variables, the condition of mass conservation of species i is given by 

vf^zfwr'' (33) 

(this equation defines the parameters /«), where 

i 

length unit also defines a dimensionless electric field via 

M = -™ = *kfZ (35) 

and a dimensionless velocity v by requiring that the relation v = /iE transforms into 

v = 6^r B v - (36) 

The diffusion coefficients A can be mapped onto length scales <2j via a Stokes formula: 
A = (37) 
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where Oj is expected to be similar to the ion radius. Nevertheless, it should be 
emphasized that the diffusion coefficient D t is a collective diffusion coefficient, not a 
tracer diffusion coefficient. With these rescalings, the Nernst-Planck equation reads 

V • (--Vc t - -^(W)c, + ~vc t ) = 0. (38) 

V a i a i J 

Finally, we introduce a dimensionless pressure via 
R 2 k B T _ 

p = W' (39) 

to re-write the Stokes equation in dimensionless form 

V •{?=(), (40) 
-Vp+ jjv 2 ^- (Vf ) J^^Cj = 0. (41) 

i 

Let us collect the final set of non-dimensionalized equations: 

V-tf=0, (42) 
- Vp + V 2 ^ - ( W) z & = 0, (43) 

i 

V • f--VQ - — Zi(V*)ci + ^c^ = 0, (44) 
V 2 ^ + ^^q = 0. (45) 

One sees that the only dimensionless parameters which remain in the equations 
are the ratios Is/cii and the charges Zj. Therefore, in order to fully characterize the 
problem, one needs to specify three parameters per ionic species (Ib/ch, Zi, and /«), plus 
the parameters which pertain to the boundary conditions: The dimensionless colloid 
radius hR, the dimensionless box size R.L (note that we assume a cubic box with periodic 
boundary conditions), and the non-dimensionalized charge density at the colloid surface. 
For the latter, we note that in conventional units the surface charge density is given by 

" = ^ (46) 

and that an electric field oriented perpendicular to the surface will jump by a value er/e. 
The jump in E is therefore given by a = ecr/(R,ekBT), i. e. 

5 =i- < 47 > 

Furthermore, (17) is straightforwardly generalized in the multi-ion case to 

$=^(^) 2 ; (48) 

this means that specification of / , Z, and RR is enough to know kL. 

We can thus summarize: In the case of zero salt and monovalent counterions, the 
reduced mobility should be a function of just the three parameters p± = RR, pi = Z, and 
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Pz — Ib/cl. This result should be contrasted with straightforward dimensional analysis, 
which was the basis of the treatment in Ref. [21]. Here one does not assume the validity 
of the electrokinetic equations, i. e. the assumption that the ionic cloud can be treated 
as a continuum field is not made. Rather, one starts from the observation that the 
problem is fully characterized by the seven parameters k B T, 77, L, Z, R, l B , and a. We 
then replace 77 by /_t (see (8)), L by R (see (29)), Z by Z (see (11)), and a by l B /a. 
This results in a new but equivalent set of parameters k B T, /i , R, Z, R, l B , and l B /a. 
Finally, we replace l B by l B /R and then R by RR to find the parameter set k B T, /i , R, 
Z, RR, l B /R and l B /a. We are thus left with seven parameters, of which k B T, /i , and 
R are needed to define the fundamental units of energy, time, and length, respectively. 
Since \i re d is a dimensionless quantity, it must be a function of the remaining four 
dimensionless parameters, which are p\ = RR, P2 = Z, p% = l B /a, and p^ = l B /R. 
Since pi, P2 and p% have also been identified on the basis of Mean Field theory (i. e. 
the electrokinetic equations), we can only conclude that any non-trivial dependence 
on pi must be the result of deviations from Mean Field theory, i. e. (most likely) ion 
correlation effects. As a matter of fact, the successful comparison between simulation 
and experimental data for fi re d that was done in Ref. [21] exclusively focused on the 
dependence on p x and p 2 - The justification for this procedure is that (i) p s is of order 
unity both in simulation and experiment, and that (ii) this implies a moderate strength 
of electrostatics. This means that ion correlation effects are expected to be weak, which 
in turn means a rather weak dependence on p 4 , and adequacy of a description in terms 
of the electrokinetic equations. 

In the case of added salt, there are further parameters which enter the problem; 
however, in the degenerate case, which was simulated in Ref. [21] and where all ion types 
have the same properties — i. e. all ions are monovalent, and have all the same mobility 
or l B /a — there are effectively only two ion types (the positive and negative ones), and 
the only additional scaling variable which enters is fo, which specifies the fraction of 
counterions relative to the salt ions. Apparently, ji re d is only weakly dependent on f 
over a wide parameter range [21]. 

In the case of finite salt concentration, we can consider the limit /o — > 0, which 
implies $ — > or L — > 00. In this case, the present formulation converges towards the 
situation studied in the "standard electrokinetic model" [6]. In the case of zero salt, it 
is not possible to perform the limit / — > within our rescaled formulation, since then 
R~ l — > 00, and this is not suitable for a length unit. However, the physics of just a free 
colloid is anyways trivial (see Sec. 1). 

Figure 1 demonstrates again the general finding that salt-free systems can, within 
reasonable approximation, be mapped onto the corresponding "salty" system with the 
same Z and RR. For a dispersion with charge Z = 20, we compare the simulation data 
for fi r ed, as a function of colloid volume fraction <3>, with the theoretical prediction that 
results from this mapping. Since it turns out that in the simulated regime of volume 
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Figure 1. Reduced electrophoretic mobility of the colloidal particle, as a function of 
the colloid volume fraction <J>. The colloid charge is Z = 20. The other parameters (in 
our simulational Lennard- Jones units) are: rj = 2.55, ksT = 1, R = 3 (friction radius 
of the "raspberry", see text), Ib = 1-3, while the coupling of the small particles to 
the lattice Boltzmann fluid is characterized by a friction constant T = 20. For further 
details on the model, see Refs. [19, 20]. The data points are simulation results. The 
solid curve is the Hiickcl formula prediction according to (49). 



fractions kR < 1, it is reasonable to assume that the Hiickel formula [1] holds: 



Here the factor (1 + RR)' 1 takes into account the reduction of the surface potential, 
within the Debye-Hiickel approximation, while Z cS was calculated via the charge 
renormalization procedure by Alexander et al. [27], based upon the Poisson-Boltzmann 
cell model. A complication arises from the fact that our "raspberry model" effectively 
defines two radii: On the one hand, the particles on the surface tether have a distance 
Ri (here: R± — 3 in our Lennard- Jones units) from the colloid center. Since the 
tethered particles are those that couple frictionally to the lattice Boltzmann fluid, this 
is the hydrodynamic radius of the sphere. On the other hand, the minimum distance 
between the ions and the colloid center is one ion diameter larger, due to the repulsive 
interaction between tether particles and ions. This defines R 2 = 4. It therefore makes 
sense to calculate the volume fraction and the RR parameter on the basis of i? 2 , and 
to also use it in the charge renormalization procedure. However, in the transformation 
from Z c fj to Z e s = Z s sIb/R, we used R±, in order to obtain the correct Stokes radius in 
the limit R — > 0. This procedure yields good agreement between simulation and theory, 
as seen from Figure 1. For the simulated $ values, Z e s varies between 16.2 and 18.4. 




(49) 



l + RR i + (3Z eff $)V2 
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4. Linear response 

In the present section, we will derive the Green-Kubo formula for the electrophoretic 
mobility, which allows us to determine /i from pure equilibrium simulations. To our 
knowledge, the formula has so far not been presented explicitly in the literature; however, 
the derivation is very straightforward within the framework of standard linear response 
theory. We follow the approach of Ref. [30], which we find particularly transparent. 
Starting point is the Hamiltonian 

H(r, t) = n (T) + n'(r, t) = w (r) - f(t)B(r), (50) 

where 7i describes the unperturbed system, and f(t) is a weak external time-dependent 
field, which couples linearly to a dynamical variable B. Y denotes the phase-space 
variable. We are interested in the dynamic linear response of a variable A. The 
time dependence of the mean value of A, A(t), must, for reasons of linearity and time 
translational invariance, have the form 

poo 

A(t) = {A)+ dr X AB(r)f(t-r), (51) 
Jo 

where (...) denotes the thermal average in the absence of perturbations. The dynamic 
susceptibility xab is defined in such a way that it is zero for negative arguments; this 
permits extension of the integration range in (51) to (—00, +00). 

For the special case that f(t) is a constant f for —00 < t < 0, and zero from then 
on, one has, for t > 0, 



A(t) = (A)+f l dr X AB(r), (52) 

or 

d 



-A{t) = -foXABit). (53) 



On the other hand, the statistical-mechanical expression for A{t) in such a "switch- 
off experiment" is 

/ drexp(-(3H + (3f B)A(t) 



A{t) fdTexpi-PHo + PfoB) ' (54) 
where (3 = l/(/c^7 1 ), A(t) denotes the time evolution of A under the influence of 7i only, 
and the Boltzmann factor describes the averaging over the initial conditions, which are 
distributed according to the perturbed Hamiltonian. Linearizing this expression with 
respect to /o for weak perturbations yields 



A(t) = f p((B(0)A(t))-(B)(A}) (55) 

or 

d 



dt A(t) = fo0(B(O)A(t)). (56) 

Comparing this with (53) yields the correlation-function expression for the dynamic 
susceptibility: 

X AB(t) = -p(B(0)A(t)) (57) 
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for t > 0. Translational invariance with respect to time implies 

= j t (B(t + r)A(t)) = (B(t + r)A(t)) + (B(t + r)i(0) , (58) 

from which one concludes, by setting r = —t, the alternative (and more useful) 
representation 



XAB(t)=P{B(0)A(t)). (59) 

Considering the case that the external perturbation is completely independent of time, 
and that A settles to a constant value, one thus finds from (51) and (59) 

A = (A) + f p dt (B(0)A(t)) . (60) 

For the problem of electrophoresis, we consider a set of particles % with charges z^e at 

— * 

positions fj, in an electric field E = Ee x . The perturbation Hamiltonian is thus given 
by 

H' = -Ee^ZiXi, (61) 

i 

i. e. f — E and B = e s ^ ji z i x i . Denoting the velocity of the ith particle in x direction 
with v ix , we thus find 

B = e^2ziV ix . (62) 

i 

On the other hand, we are interested in the velocity response of one particular particle 
(say, % = 0), i. e. 

A = v 0x , (63) 

with (A) = 0. This yields directly the desired Green-Kubo formula for the electro- 
phoretic mobility 

1 e v ^ f°° 

^=Sk^T^ Zt J * <*(°) , (64) 

where we have averaged over the three spatial directions. It should be noticed that 
the formula involves a mixed correlation between the test particle and all charges, in 
contrast to the tracer diffusion coefficient, which contains only the autocorrelation of 
the test particle, and the electric conductivity, which involves the autocorrelation of the 
collective current. 

As an example, we present Figure 2, where the reduced mobility for a salt-free 
system is plotted as a function of colloidal charge. Comparison with the Green- 
Kubo integral makes it possible to check whether data obtained under non-equilibrium 
conditions are still within the linear regime or not. One sees that the mobility first 
increases with the charge (as one would expect from the physics of the free colloid), but 
then saturates at a finite value, as a result of condensation of more and more counterions. 
The nonlinear effects observed for stronger electric fields are mainly due to charge-cloud 
stripping [20], which increases the effective charge and thus the mobility. 
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Figure 2. Reduced electrophoretic mobility of the colloidal particle, as a function of its 
charge Z. No salt is added, and, apart from the central colloid, the system comprises 
Z monovalent counterions. The linear box size is L = 30. The other simulation 
parameters are the same as those given in Figure 1. The mobility was here defined 
just as the ratio between drift velocity and electric field. For strong driving field E, 
one observes nonlinear effects, while the results for weak driving agree favorably with 
the results of Green-Kubo integration (GKI). Note that the driving field is here given 
in the Lennard- Jones unit system of the simulation. According to (35), constant E 
does not imply constant E, since R varies with Z . 



5. Finite element calculations 

As a complementary approach to the hybrid MD / LB simulations, we have also 
done some calculations where we solved the electrokinetic equations directly, using 
a commercial finite-element software package (COMSOL 3.3). For highly charged 
systems, where rather fine grids are necessary, this does not work particularly well, 
since quite generally the software tends to need excessive amounts of memory. We used 
the same geometry as in the simulations, with the colloidal sphere centered in the cubic 
box, but confined, for simplicity, the computational domain to just the space outside 
the colloidal sphere. This is not entirely correct, since, in reality, the electric field also 
exists inside the sphere, where it takes a non-trivial configuration, as a result of the 
external driving field oriented in x direction, the deformed charge cloud, and the cubic 
anisotropy. However, if we assume that we can neglect the latter, and consider the limit 
of infinitesimal driving, we get an electric field at the colloid surface whose orientation 
is strictly radial, and whose value is given by Gauss' law. This corresponds to the 
specification of Neumann boundary conditions for the normal component of the electric 
field. On the surface of the box, we specified Neumann boundary conditions as well, 
where the normal component was set to zero in the planes perpendicular to y and z, 
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Figure 3. Reduced mobility as a function of electric field E, for parameters chosen 
in accordance with those of Figure 2, for charge Z — 60. Note that the electric field 
is given here in reduced units (see (35)), i. e. a value of E = 1 in the present plot 
corresponds to eE w 0.2 in the Lennard- Jones units of Figure 2. 



while on the planes perpendicular to x it was set equal to the driving field. Concentration 
and flow field were subjected to periodic boundary conditions. The pressure and the 
electrostatic potential were set to zero at some arbitrary point in the domain, in order 
to lift the degeneracy of shifting these functions by an arbitrary amount. The Nernst- 
Planck equation was augmented by a zero-flux condition at the colloid surface, and an 
integral constraint in order to guarantee charge neutrality (such integral constraints turn 
out to be computationally particularly cumbersome). The flow velocity at the colloid 
surface was set to zero, and the particle velocity was finally determined by transforming 
back into the system's center-of-mass reference frame. Given the inaccuracies of the 
boundary conditions, these results should not be viewed as a stringent test of the validity 
of the Mean Field picture for the simulated system. Nevertheless, the results agree 
reasonably: Figure 3 shows the reduced mobility as a function of the external driving 
field, for a situation which corresponds to the parameters of Figure 2 at charge Z = 60. 
In the future, we hope to be able to calculate reduced mobilities in the low-salt limit 
more easily by self-written software; efforts to develop such a program are currently 
under way. 

6. Summary 

In this paper we developed a theoretical basis for the scaling analysis of the colloidal 
electrophoresis problem in the case of finite colloidal concentrations. The rescaling 
procedure and characterization of the dispersion in terms of effective dimensionless 
parameters, i. e. the reduced colloid charge, and the ratio of screening length and 
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size, allows one to map the numerical results obtained for a single colloid onto 
experimental data for finite colloidal volume fractions and no added salt. At least 
for a certain parameter regime, we can also map the salt-free case onto a corresponding 
system containing additional salt. Moreover, we presented a numerically convenient 
method of measuring the colloidal electrophoretic mobility based on the Green-Kubo 
analysis of the equilibrium fluctuations of the charge motions. This allows for pure 
equilibrium simulations and ensures that one always measures the mobility and ion 
distributions in the linear response regime. Finally, we gave an example of using 
a finite-element commercial software package for solving the electrokinetic equations 
numerically, yielding reasonable agreement with the simulations, and suggesting at least 
consistency of the Mean Field picture with our simulations. 
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